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Abstract 

Recent density functional theory (DFT) calculations by Forst et al^ have predicted that va- 
cancies in both low and high carbon steels have a carbon dimer bound to them. This is likely to 
change the thinking of metallurgists in the kinetics of the development of microstructures. While 
the notion of a C2 molecule bound to a vacancy in Fe will potentially assume a central importance 
in the atomistic modeling of steels, neither a recent tight binding (TB) model nor existing classical 
interatomic potentials can account for it. Here we present a new TB model for C in Fe, based on 
our earlier work for H in Fe, which correctly predicts the structure and energetics of the C2 dimer 
at a vacancy in Fe. Moreover the model is capable of dealing with both concentrated and dilute 
limits of carbon in both a-Fe and 7-Fe as comparisons with DFT show. We use both DFT and 
TB to make a detailed analysis of the dimer and to come to an understanding as to what governs 
the choice of its curious orientation within the vacancy. 

PACS numbers: 71.20.Be 75.50.Bb 73.20.Hb 68.43.Fg 



1 



I. INTRODUCTION 



A great deal of progress over many years has been made in understanding the physics 
of the bonding and electronic structure in pure magnetic iron^ and this has been used with 
great effect to advance the generation of schemes from density functional theory in the gra- 
dient corrected^ local spin density approximation^ (LSDA-GGA) through the tight binding 
approximation 6 -^ to classical interatomic potentials to be used for atomistic simulations by 
the materials science community.™ Of course the metallurgist is rarely interested in pure iron 
and so the challenge to the physicist has been to extend the theory to include interstitial 
carbon, which is the defining element whose presence distinguishes steel from iron. It is 
only in the last decade that real progress has been made; first with some very extensive 
LSDA-GGA calculations,-*^— second with the generation of (admittedly very complicated) 
classical interatomic potentials based in the embedded atom method,— 1 ^ and third, the 
subject of this paper, by some recent semi empirical quantum mechanical schemes based 
in the tight binding approximation. A particularly significant advance was made recently 
by Hatcher, Madsen and Drauta^ who constructed a very simple orthogonal tight binding 
model for carbon and iron using a minimal basis of C-p and Fe-d orbitals and a local charge 
neutrality condition. This model is a natural basis for a bond order potential, 16 but we 
argue here that this basis may be too small to capture some of the physics of carbon in iron. 
Instead we introduce a new model based in our earlier work on H in Fei£ employing a larger, 
non orthogonal basis of C-p, Fe-d, and C and Fe-s orbitals and treating charge transfer self 
consistently via an adjustable "Hubbard-t/" parameter.— The structure of the paper is as 
follows. In section [Til we describe our new model for carbon in iron and in section II III we 
demonstrate its predictive power in both the concentrated (iron carbide) and dilute impurity 
limits. In lHICI we focus on the carbon dimer bound to a vacancy, taking in view the recent 
startling prediction from LSDA-GGA that this is a predominant point defect in steel.- Our 
discussion and conclusions are to be found in sections [TV] and [V] 

II. DESCRIPTION OF THE NEW MODEL 

We take the same approach as in our earlier work on H in Fe¥- which is similar to that of 
Hatcher et a/.— on C in Fe, namely to proceed from a given model for pure Fe and generate a 
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further parameterization for the interstitial element. In contrast to Hatcher et al— we do not 
use a direct projection of the LSDA-GGA Hamiltonian onto a tight binding basis,— >2S instead 
we employ a genetic algorithm^ to fit the parameters to a small set of LSDA-GGA target 
data which enter an objective function, which is minimized. As a consequence of employing 
the simplest tight binding scheme, namely an orthogonal basis of only <i-orbitals on the 
Fe atoms and p-orbitals on the C atoms, the model of Hatcher et al.— differs significantly 
from ours. One difference results from their underlying model for pure Fe which includes 
an attractive bonding term in the total energy which is environment dependent and which 
accounts for a significant fraction of the total energy— This was intended as a surrogate for 
the missing s-electrons, but the fact that this term is large and negative is surprising as one 
expects the s-band to exert a positive pressure.— ~— Another difference is that in the minimal 
basis having only p-orbitals on C atoms the limit of pure carbon can only be approximately 
rendered since it is the sp-hybridisation in carbon that leads to the rich variety of single, 
double and triple bonds and the competition between sp 2 -bonded graphite and sp 3 -bonded 
diamond. We will argue below that carbon sp-hybridization plays a key role in the structural 
stability of iron carbides and also in controlling the configuration of the C2 dimer at an Fe 
vacancy. Therefore in the current work we employ a larger basis, namely s- and d-orbitals on 
Fe and s- and p-orbitals on C atoms from which we suppose that at the expense of greater 
computational cost we have a physically better motivated model. Moreover we use a non 
orthogonal basis, and as we argued earlier— we believe that this allows a more natural way 
to include environment dependence in the bond energy. 

The tight binding model that we present here is identical in its mathematical form to 
those we developed earlier.— The functional forms of the bond and overlap integrals are 

h{r) = h e- qr s{r) = s e- qr 

and we tabulate all parameter values ho and so in tables [J (for Fe-Fe terms) and |V] (for Fe-C 
interactions). The Fe-Fe and Fe-C pair potentials are 

0(r) = B x e~ pir - B 2 e~ P2r 

noting the sign, so that in tables HI and IVl parameter values for Bi and B 2 are positive. For 
Fe-H the pair potential is 




TABLE I. Intersite bond integral and pair potential parameters for the Fe-Fe terms in our tight 
binding model. All quantities are in Rydberg atomic units (a.u.) except for the cut off radii, t\ and 
r c which are in units of the equilibrium bcc lattice constant aQ CC =2.87A. We include models using 
both volume-scaled and fixed cut offs (see the text). These differ only in their pair potentials and 
are indicated as "sc" and "fx" respectively in the last four rows. Properties of pure Fe resulting 
from these four models are displayed in table HU Parameters for C and H in Fe in table [V] and used 
in the remainder of this paper are associated with the fixed cut off sd- model (scZ-fx). We note that 
these parameters differ from those published earlier; 17 first by correcting misprints in the decay 
constant q in the ssa and sda terms, second because we have moved the cut off r\ from 10% larger 
to 10% smaller than ajf c . Third we now prefer fixed multiplicative to scaled augmentative cut offs. 
See also the text. These differences are then reflected in slightly different calculated properties in 
table H 



ssa 


sda 


dda 


ddir 


dd5 


ho s 


ho s 


h 


h 


ho 



-0.35 0.45 -0.14067 0.5 -2.4383 1.9972 -0.90724 



_q = 0.3 q = 0.9 

n = l.i n = 0.9 

r c = 2.0 r c = 1.4 





Bi 


Pi 


B 2 


P2 


n 


r c 


sd-£x 


698.67 


1.52 


517.467 


1.4576 


0.9 


1.4 


sd-sc 


665.60 


1.40843 


536.800 


1.36297 


0.9 


1.4 


d-fx 


683.1 


1.5376 


459.5 


1.4544 


0.9 


1.4 


d-sc 


682.8 


1.5165 


466.8 


1.4350 


0.9 


1.4 
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TABLE II. Properties of pure a-Fe. Target elastic constants are taken from experimental low 
temperature data (see ref theoretical hcp-bcc energy difference from our own LSDA-GGA 

calculations; experimental vacancy formation and migration energies are from Seeger,— while re- 
maining LSDA-GGA data are from Domain and Bequart,— and Kabir et alM- We show four models 
in data columns 1-4: canonical (d) with scaled (sc) and fixed (fx) cut offs and non orthogonal (sd) 
with scaled and fixed cut offs respectively. 





d-sc 


d-fx 


sd-sc 


sd-fx 


target 




K (Gpa) 


161 


174 


185 


192 


168 


(expt.) 


C (GPa) 


50 


50 


55 


45 


53 


(expt.) 


C44 (Gpa) 


118 


117 


106 


100 


122 


(expt.) 


A^kp-bcc (mRy) 


8 


6 


6 


6 


15 


(LSDA-GGA) 


<c. (eV) 


2.0 




1.6 




1.61-1.75 
2.0 


(expt.) 
(LSDA-GGA) 


HfL ( eV ) 


1.16 




0.81 




1.12-1.34 
0.65-0.75 


(expt.) 
(LSDA-GGA) 



We take a rather sophisticated approach to cutting off the spacial dependence of these 
interactions. We require proper energy conservation in molecular dynamics and cannot allow 
discontinuities in second derivatives of bond integral or pair potential functions. Previously^ 
we implemented the cut off by augmenting (that is, replacing) the function with a polynomial 
of degree five within rq < r < r c whose coefficients are chosen so as to match the function 
continuously and differentiably to its value at rq and to zero at r c . We chose rq = Llajj 00 and 
r c = 1.4aQ CC , where Oq cc = 2.87 A is the lattice constant of a-Fe, so that functions are cut 
off to zero between second and third neighbors of the bcc lattice. Subsequent improvements 
were introduced after making two observations.— (?) A smoother effect can be achieved 
using a multiplicative cut off; that is, to multiply the function by a polynomial of degree 
five whose value is one at rq and zero at r c and whose coefficients again ensure that the 
function is everywhere continuous up to the second derivative, (ii) Because the multiplicative 
cut off "inherits" the shape of the function near rq better than the augmentative cut off, 
we found that we could move rq back to 0.9aQ CC and achieve a smoother function overall. 
A second difference compared to our earlier work— is that there we employed a volume 
dependent cutoff, whereas now we prefer to use a cut off that is fixed. Because of this 
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TABLE III. On-site Hamiltonian matrix elements for tight binding models of Fe, H and C. U is 
the Hubbard-?/ parameter and J is the Stoner parameter.— All quantities are in Rydberg atomic 
units (a.u.); Nd is the number of (^-electrons, which is an adjustable parameter in the canonical 
model. 

e s -ef e p -ef N d U J 
7 0.05 

0.15 1.0 0.055 

0.468 0.083 1.238 

0.085 1.2 



Fe-d 
Fe-sd 
C 
H 



small modification it is necessary to obtain slightly amended pair potential parameters. We 
show these in table HJ and in table [III some predicted properties of pure Fe using both the 
canonical <i-band model for Fe and the non orthogonal sci-model. Our canonical model can 
be read from table [J simply by ignoring those parameters that don't enter the Hamiltonian. 
This model therefore differs from the canonical model that we published earlier.— In fact 
it is worth pointing out that both these canonical d-band models reproduce the vacancy 
formation and migration energies better than our non orthogonal model; although its -ffyac 
is outside the experimental range it is in better agreement with published LSDA-GGA data. 
In our opinion the simplest canonical model is very appropriate for pure Fe and we do not see 
the need for the additional, attractive "embedding potential" introduced by Madsen et al— 
The inclusion of the s-band and non orthogonality in Fe is only necessary once hydrogen or 
first row elements are included. The reason for this is that the valence s-band from these 
elements lies typically below the Fe 3(i-bands; orthogonality constraints in the concentrated 
limit then push the Fe 4s-band to above the Fermi level. In the dilute limit the electronic 
structure has to differentiate between regions close to an impurity and those far from it 
where the iron 4s local density of states returns to its position in pure Fe below the Fermi 
level. To account properly for this effect the impurity and Fe s-bands cannot be neglected. 
Of course in a minimal pd basis for Fe-C tight binding models or bond order potentials both 
s-bands are neglected which is internally consistent, but these models cannot account for, 
say, the carbon sp-hybridisation. 
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FIG. 1. (color online) Structural energy- volume curves for the four iron monocarbide phases, FeC, 
that were used in the fitting of the tight binding model. These show the heat of formation as 
a function of the atomic volume per Fe atom in both bcc (a) and fee (7) Fe each containing C 
atoms in tetrahedral (TET) or octahedral (OCT) interstices. Note that the atomic volume, f^o, 
of pure a-Fe is 79.765 a.u. In this, and subsequent figures, GGA denotes the generalized gradient 
approximation to the LSDA. 
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TABLE IV. Properties of iron monocarbides. These are the atomic volume of Fe, f2p e5 in FeC 
relative to $7o = ( a o cc ) 3 /2 in the four monocarbides obtained from bcc a-Fe and fee 7-Fe in which 
C is in either tetrahedral or octahedral interstices. AE is the energy difference per formula unit in 
Ry compared to the fee octahedral compound. Targets are taken from our LSDA-GGA calculations. 





a-tet 




Q-OCt 




7-tet 




7-oct 




^Fe/^0 


AE 


^Fe/^0 


AE 


^Fe/^0 


AE 


^Fc/^0 


TB 


1.591 


0.030 


1.832 


0.202 


1.690 


0.005 


1.327 


target 


1.549 


0.020 


1.788 


0.147 


1.613 


-0.010 


1.339 
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A. The on-site carbon and Fe — C parameters 

Our approach to finding carbon on-site energy parameters, Fe-C Hamiltonian matrix 
elements and pair potential parameters is to fit these to just seven target data (a refinement 
was done later to improve the model in the dilute limit, see section UlI Bl) . These data are 
taken from LSDA-GGA calculations and illustrated in figure [I] which shows the heat of 
formation of four compounds having the stoichiometry FeC. These are either bcc a-Fe or 
fee 7-Fe with C interstitials in tetrahedral or octahedral sites. In the 7-Fe case these are 
identical to the zincblende and rocksalt crystal structures. We fitted our parameters to 
the four equilibrium volumes and three energy differences. The outcomes of the fitting are 
shown to the left of figure [1] and in table IIVI The resulting parameter values are tabulated 
in tables (TJ 1111} and |VJ For comparison and for completeness we also show parameters of 
our earlier model for hydrogen in Fe- 
lt is important to make some comments about the energy-volume curves in figure HJ also 
in relation to the equivalent data for the hydrogen interstitial.— First, in the bcc structure 
both C and H prefer the tetrahedral site in the concentrated limit of FeC and FeH, and 
this remains the preferred site for H into the dilute limit. In contrast carbon occupies the 
octahedral sites in both ferritic and austenitic steel. In a-Fe, this is achieved at the expense 
of a local tetragonal distortion of the lattice so as to drive apart the two apical Fe atoms 
in the irregular octahedron of the underlying bcc lattice. This is only possible if the C is 
sufficiently dilute, certainly more dilute than the stoichiometry Fe4C, as we will see below, 
and in fact the crossover is around FeigC— Second, in the fee structure it is certainly striking 
that according to LSDA-GGA FeC adopts the zincblende structure rather than the rocksalt 
structure, albeit at an expanded volume, as the tetrahedral interstice is much smaller than 
the octahedral. This is contrary to the behavior of hydrogen, even though its atomic radius 
is evidently smaller. Again there is a crossover towards the dilute limit where C prefers the 
octahedral site in 7-Fe.— We expect that the competition between the two sites in FeC is 
driven by the sp-hybridization which will be maximal in the four fold coordinated tetrahedral 
site, whereas the six fold octahedral site offers a bonding environment favorable to the 90° 
bond angles of unhybridised p-orbitals. Therefore it is surprising that the p<i-basis model 15 
reproduces this result correctly. We expect that this arises from the freedom of employing 
long ranged C-C interactions in that model.— Conversely we take the canonical point of view 
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TABLE V. Intersite bond integrals and pair potential parameters for the Fe-C and Fe-H terms. 
All quantities are in Rydberg atomic units (a.u.) except for the cut off radii, n and r c which are 
in units of the equilibrium bcc lattice constant <2q cc = 2.87A. This corrects two misprints regarding 
the ssa and sda bond integrals for Fe-H in table V of ref [171 ]. 
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Fe-C 


-1.7712 


0.38434 


3.9546 -0.59202 


-0.17549 


0.10283 


-1.2300 0.32895 


0.88500 -0.37025 




0.56548 


0.30106 


0.76024 0.39114 


0.30249 


0.34080 


0.64362 0.30636 


0.66529 0.45518 




n 598 


n 611 






595 






1 790 


1 644 






1 674 




Fe-H 


-1.0935 


0.26587 




-0.40748 


0.21988 








0.77628 


0.28633 




0.45450 


0.47301 








0.8 


0.8 




0.8 


0.8 








2 


2 




2 


2 












B 1 p 1 


B 2 


P2 


r\ r c 








Fe-C 


771.190 2.3962 


19.325 


1.5555 


0.50071 1.5070 








Fe-H 


299.563 2.69225 






0.75 0.95 





that C-C interactions should not extend beyond the first neighbor distance in diamond, 28 
as it is known that longer range terms do not improve tight binding models for diamond 
structure sp-bonded elements.— 1 ^ 



B. The C — C parameterization 

Taking the view that carbon-carbon interactions are to be curtailed beyond the usual 
definition of the chemical bond lengths of 1.2-1.5 A, none of the tests that we will apply 
in section II I II will require us to specify the C-C bond integrals or pair potential. There is 
however one notable exception which we will be discussing in greater detail below. This is the 
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observation from LSDA-GGA calculations^ that two carbon atoms bound to a monovacancy 
in Fe will form a "dimer molecule" whose bond length is 1.44 A. In our model we have the 
freedom to choose our C— C interactions at will since they do not affect any of the results 
in the concentrated limit. It would be desirable if an existing model for diamond could be 
adopted without modification and we have used a tight binding Hamiltonian for diamond 
from Harrison^ with parameters adapted by Xu et al. 28 (Our model is essentially that of 
Harrison in terms of the scaling of the bond integrals and pair potential with bond length. 
We take over the bond integrals at the equilibrium volume in diamond from Xu et al. but 
we do not adopt their scaling.) In this way for the C-C parameters we use a simple power 
law model, namely 

h(r) = h Q r~ 2 , 0(r) = r~ 4 

with (in Rydberg a.u.) 

h s Q sa = -3.734 h s pa = 3.510 h p pa = 4.107 ^ = -1.157 

B\ = 50 a.u. leads to the correct lattice constant and bulk modulus in diamond carbon. 
Unfortunately that choice of B\ does not quite reproduce the energy and C-C bond length 
of the C2 molecule bound to the Fe monovacancy. Therefore we employ B\ = 43 a.u. which 
leads to about an 8% error in the diamond lattice constant. We have used the same value of 
Bi to calculate the total energy of diamond which is the quantity we have used in figures dHU 
to determine the heat of formation of Fe-C compounds from elemental a-Fe and diamond C. 
This has the consequence that the tight binding theory consistently overestimates —AHf 
by 0.16 Ry (see for example figured]); if we use the value B\ = 50 a.u. the agreement with 
the LSDA-GGA is rather better. To keep the C-C interactions to within the first neighbors 
as expected in diamond and in hydrocarbons, we apply a fixed multiplicative cut-off at 
n = 0.6<2q cc and r c = cIq cc . 

III. PREDICTIONS OF THE NEW MODEL 

In this section we examine to what extent our model reproduces some previously published 
or our calculated LSDA-GGA data. 
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FIG. 2. (color online) Energy volume curves for compounds with stoichiometry Fe2C, comparing 
the predictions of our model with results of our LSDA-GGA calculations. These compounds have 
either bcc (a) or fee (7) iron lattices with C placed at interstitial positions in tetrahedral (t), 
octahedral (o), in the bcc case the saddle point (s) along the (110) direction and, in the fee case, 
the saddle points (s) along the (111) direction and (d) along the (110) direction. This latter site 
is midway between two nearest neighbor Fe atoms and so is expected to have a high energy; on 
the other hand as is known from LSDA-GGA calculations 27 and as our model also predicts, this 
site is along the diffusion path in 7-Fe. The alternative path for diffusion (adopted by H) via an 
intermediate tetrahedral site has higher energy. This is discussed in section IIIIB1 




£2 Fe / a.u. £2 Fe / a.u. 

A. The concentrated limit 

We first compare tight binding with LSDA-GGA for a range of mostly fictitious Fe-C 
compounds with stoichiometries Fe2C, FesC, and Fe4C in figures [2HU The LSDA-GGA 
calculations were done by means of the mixed-basis pseudopotential (MBPP) method.— 
The PBE-GGA exchange and correlation functional,- optimally smooth norm conserving 
pseudopotentials^ for Fe and C, k-points which are equivalent to 8 x 8 x 8 Chadi-Cohen 
meshes for cubic structures, and a Gaussian broadening by 0.2 eV were employed. The mixed 
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FIG. 3. (color online) Energy volume curves for compounds with stoichiometry Fe3C, comparing 
the predictions of our model with results of our LSDA-GGA calculations. Here we show both 
interstitial, and substitutional putative phases all of which have large positive heats of formation 
and are hence predicted not to exist. The first four in the column of labels represent bcc and fee 
supercells containing one vacancy and a C atom at either a neighboring tetrahedral (t) or octahedral 
(o) site. The next four are substitutional phases labeled using their Strukturbericht designations. 
Of greatest interest are the remaining three FesC compounds: "WC" labels a hypothetical high 
energy structure similar to a simple hexagonal tungsten carbide like phase (see the text); the only 
phases predicted to exist thermodynamically are the e carbide and the 9 carbide, or cementite 
phase. It is interesting to note that the LSDA-GGA predicts these both to have a positive heat of 
formation, which requires further investigation, since they are both known to exist ubiquitously in 
the microstructures of steels. 




40 60 80 100 120 60 80 100 120 

Q. Fe / a.u. Q. Fe I a.u. 



12 



FIG. 4. (color online) Energy volume curves for compounds with stoichiometry Fe4C, comparing 
the predictions of our tight binding model with our LSDA-GGA calculations. Here we show a 
variety of supercells in bcc (a) and fee (7) crystal structures, each in face centered cubic (fee), 
body centered tetragonal (bet) and simple cubic (sc) supercell settings. In each of these cases, a C 
atom is placed at one of the interstitial sites, using the same labeling as in figure [2j As expected 
from the results of figure [T] and seen also in the data in figure [2] only the 7-Fe with octahedral 
carbon have reasonably low heats of formation. 



\\ 1 I 1 I 1 I 1 I §~E 

E 



fit 0.3 



SB 



1 r~p I ' I ' I 1 1 r- 




100 110 120 
Q. Pc I a.u. 



100 110 120 
£2 Fe / a.u. 



basis consisted of plane waves up to a maximum kinetic energy of 340 eV and atom-centered 
basis functions with (i-symmetry for Fe atoms and with p-symmetry for C atoms which 
are confined to spheres with radii of 1.1 9 A and 0.66A respectively. The broad agreement 
between tight binding and LSDA-GGA is excellent and indeed in almost all instances the 
ordering in energy of the phases is correctly reproduced. 

We will focus most closely on the stoichiometry FesC, figure [3j which is the most signif- 
icant composition in materials science due to the ubiquitous occurence of cementite in the 
microstructures of steels. It is also notable that a significant weakness in the tight binding 
model emerges here when trying to describe the hypothetical substitutional phases D0 3 , L6 , 
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TABLE VI. Calculated lattice parameters of hexagonal e FesC and orthorhombic Fe%C (cemen- 
tite), compared to LSDA-GGA calculations and experimental measurements reported by Jang et 
al l 34 ' 35 Lattice parameters a, b and c are in A. The last column shows the predicted equilibrium 
volume compared to experiment. 







a 


b 


c 


b/a 


c/a 


v/v cxp 


0-Fe 3 C 


TB 


4.95 


6.79 


4.42 


1.37 


0.89 


0.96 




LSDA-GGA 


5.13 


6.65 


4.46 


1.30 


0.86 


0.98 




exp. 


5.09 


6.74 


4.52 


1.32 


0.89 




e-Fe 3 C 


TB 


4.63 




8.64 




1.87 


0.93 




LSDA-GGA 


4.74 




8.63 




1.82 


0.98 




exp. 


4.77 




8.71 




1.83 





LI2, and D022- This is not so surprising since this bonding environment is very different 
from the interstitial phases and fortunately the substitutional phases are not of particularly 
great interest. It is on the other hand very gratifying that the tight binding model repro- 
duces with great fidelity the e and 9 iron carbide phases. At the same time the ordering of 
the unfavorable simple hexagonal tungsten carbide like structure is very well rendered; this 
hypothetical structure is obtained from the hexagonal close packed e-FesC by rotating alter- 
nate Fe layers about the c-axis by 60° and increasing the axial ratio by a/3/2. To emphasise 
the suitability of the tight binding approximation in the modeling of steel microstructures, 
we show in table lVIl a detailed comparison of the calculated crystal lattice parameters of the 
important phases, e and 9 iron carbide, with experimental data. 

B. The dilute limit 

Of equal or even greater interest is the behavior of carbon in the dilute limit. In the 
case of hydrogen in Fe we could claim a success in that a model fitted in the most con- 
centrated stoichiometry transfers very well into the dilute limit. 1 - Carbon in Fe has been 
more problematic and subsequent to the fitting described in section Hi A I it was necessary to 
make further genetic optimizations of the Fe-C parameters in order to render correctly the 
migration energy and the binding of C to a monovacancy — the quantities H^} a and Eb(1) 
in the first two data columns of table IVII1 
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TABLE VII. Properties of C in Fe in the dilute limit. H^ a is the migration energy of the C 
atom, equal to the energy difference between C in tetrahedral and octahedral interstices in bcc 
a-Fe. Eb(1) and -Ee(2k) are the binding energies of one and two C atoms to a vacancy (see 



ref 



111 ], table 5). The final two columns are data for C in fee 7-Fe and show the migration energies 



H^"' (tet) and (d) of C between two octahedral sites via a tetrahedral site and the "d-sadd 



site respectively. All energies are in eV. The first line shows LSDA-GGA data taken from refs 9j 
H and fl. The second line shows results from our TB model. 





H c 


E B {1) 


E B {2k) 


F^(tet) 




LSDA-GGA 


0.87 


0.47 


1.50 


1.48 


1.00 


TB 


0.81 


0.35 


1.55 


2.11 


0.63 



TABLE VIII. Comparative energies and bond lengths of four possible configurations of two carbon 
atoms bound to one vacancy in bcc a-Fe. We use the designations of refs A and Q. The 
structures relaxed using TB are illustrated in figure [5] The first and second rows show results 
from LSDA-GGA calculations and the third and fourth those from the TB model. The latter 
correctly predicts the ordering in energy (these are shown relative to the "k" ground state energy) 
and the bond length. The "(110)" configuration has the dimer centered at the vacant Fe site and 
orientated along a (110) direction. Energy differences AE are in eV and bond lengths d in A. The 
numbers in parentheses are energy differences calculated using the Fe atom positions of the relaxed 
equivalent structure, with both C atoms removed; these numbers allow a comparison of the host 
lattice distortions accompanying the formation of the dimer-vacancy complex. 

"j" "(100)" "k" "(110)" 

AE d AE d AE d AE d~ 

LSDA-GGA 0.37 2.57 0.11 1.46 1.43 0.06 1.43 

(0.07) (-0.05) (0) (0.04) 

TB 1.10 2.70 0.13 1.46 1.44 1.23 1.40 

(0.49) (0.16) (0) (0.66) 
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Our main results are presented in tables IVHI and IVIIIl We have constructed cubic 4 x 
4x4 supercells for a-Fe and 7-Fe in order to study, in particular, the energetics of the 
monovacancy in Fe and its binding to carbon interstitials. As in the case of hydrogen, the 
impurity does not occupy a vacant Fe lattice site, as is clear from figure [3] which shows a 
large positive heat of formation for the four substitutional phases considered. Instead (again, 
as does hydrogen) carbon occupies a position close to its preferred interstitial site, in this 
case the octahedral interstice, at one of the cube faces bounding the vacancy. We follow 
the definitions employed by Becquart et al— such that the binding energy of one or more 
interstitials to a vacancy is the difference in energy between that number of interstitials and 
the vacancy occupying separate, non interacting sites, and the interstitials bound to the 
vacancy. In this way, we have, from calculations based on a 128-atom supercell of pure Fe, 

E B {1) = -(£(Fe 127 C) + £(Fe 128 )) + (E(Fe 127 ) + E(Fe 12S C)) 

and 

E B (2) = -(£(Fe 127 C 2 ) + 2£(Fe 128 )) + (£(Fe 127 ) + 2£(Fe 128 C)) 

where the signs are employed such that a positive binding energy implies a preference for 
the two C atoms to bind at a vacancy compared to the vacancy and two C interstitials 
being widely separated. The total energies E involved are calculated by relaxing supercells 
containing the numbers of atoms indicated in parentheses. 

We therefore show in table rvTIl -EWl). the binding energy of a single C atom to a mono- 
vacancy, and Eb (2k) (using the designations of Becquart et al.—) the binding energy of two 
C atoms to a vacancy. 

We have also calculated migration energies of carbon in a-Fe and 7-Fe using static re- 
laxations and also the nudged elastic band method. 3 - TB describes these correctly in both 
phases of Fe as seen in the data columns 1, 4, and 5 in table IVHI In particular our TB 
model confirms the LSDA-GGA result that the diffusion path of C in 7-Fe is not as one 
might suppose mediated by a "double" hop via a neighboring tetrahedral site as for H in 
7-Fe, but the carbon atom actually takes a direct route forcing itself through the bond center 
of two nearest neighbor Fe atoms at the (110) (d) saddle point. This is surprising in view of 
the high energy of the 7-d crystal structures in figures [3] and HI However this tight binding 
prediction agrees with LSDA-GGA results.— 
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C. The carbon dimer at the vacancy 



Several authors have made LSDA-GGA calculations for a carbon dimer bound to a va- 
cancy in a-Fe and described a number of possible atomic structures.— liiiM We consider four 
here. If the dimer is orientated along a (100) direction, then if the carbon atoms remain 
close to their original octahedral sites at opposite faces of the cube bounding the vacancy, 
this configuration is designated "j" by Becquart et alM or "00" by Lau et al— This is a 
local minimum in the energy and in our TB model the two carbon atoms are outside their 
range of interaction in the Hamiltonian (see section HlBj) . The energy is lowered if the two 
carbon atoms approach each other along the (100) direction and form a dimer bond, having 
a bond length of 1.46A which is very close to that in diamond and the C-C single bond in 
molecules. This configuration is denoted "(100)" by Lau et al.— but was not considered in 
earlier work.— »ii This is not yet the global minimum energy for the dimer which is achieved 
by orientating the dimer along a (Oil) direction with the two C atoms close to octahedral 
positions, a situation denoted "k" by Becquart et al— or "AO" by Lau et al— The latter 
authors describe a further configuration, "(Oil)" in which the dimer bond is centered at the 
vacant site and orientated parallel to "k". This is of slightly higher energy than "k". The 
four configurations are are illustrated in figure In table IVIIII we compare predictions of 
the TB model with our own LSDA-GGA results.— The only serious discrepancy is that the 
TB model overestimates the energies of the "j" and "(Oil)" configurations with respect to 
LSDA-GGA. 



IV. DISCUSSION 

Forst et al— made a very thorough study using LSDA-GGA of point defect complex 
energetics and found the remarkable result that under the conditions normally encountered 
in a steel, effectively all Fe vacancies have a C2 dimer bound to them as illustrated in figure O 
Moreover, by just a few hundredths of an eV, the dimer prefers to be orientated along a 
(110) direction. It would be of great interest to determine whether this phenomenon can 
be confirmed experimentally, possibly by internal friction measurements. It is notable that 
a similar prediction was made using DFT concerning dimerization of boron in copper, a 
prediction that is consistent with thermodynamic assessment, 41 
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FIG. 5. (color online) Atomic structures of four configurations for a carbon dimer bound to a mono- 



vacency in a-Fe. Structures shown are (a) "j" [111], or "00" 



11 ] or "AO" []j}], the global minimum for this configuration. 14 The relaxed structures displayed 



14]; (b) "(100)" 



14j];(c)"(110)"[l4j];(d) 



here are obtained with our TB model; Fe-C bond lengths shown are 3.73A in "(100)" and 3.65A 
in "k" 




18 



One can interpret the competition between the four configurations illustrated in figure [5] 
in terms of certain notional contributions to the total energy. These are (i) the formation 
of a C— C covalent bond, (ii) the coordination of the carbon atoms to neighboring Fe atoms, 
and (iii) the accompanying distortion of the host Fe lattice containing a vacancy, (i) It is 
surprising that a C-C bond having the same length as in pure carbon is formed in view 
of the metallic electron gas destroying the single bond order; the bond length inside the 
metal is the same as in the pure carbon or hydrocarbon, but the bond energy is about 
ten times smaller, (ii) The Fe-C coordination goes a long way to explain the stability of 
the most stable configuration "k" in which the carbon dimer makes three bonds of equal 
length (3.65A) to Fe atoms, thus forming an "ethane" molecule in which the Fe atoms take 
the part of H atoms (see figure E(d)). With the dimer orientated along (100) each carbon 
atom makes four bonds (3.73A long) to neighboring Fe atoms (figure EJ^b)). We take it that 
this is less favorable owing to carbon preferring a four fold coordination. Configurations 
"j" and "(110)" both display a planar configuration of Fe-C bonds; in "j" the carbon is 
bonded to four Fe atoms in the plane of a cube face of the bcc lattice, in "(110)" each 
carbon is bonded to two Fe atoms, (iii) In table [yTTTI we give in parentheses values of the 
calculated total energies of the four configurations, having removed the two carbon atoms and 
leaving the Fe atoms in their positions. This is intended to examine the elastic distortion 
energy accompanying the introduction of the dimer into the vacancy. In the LSDA-GGA 
these distortion energies are rather small and in fact "(100)" has a slightly lower value than 
"k". However the preferred four fold coordination of the carbon atoms in "k" is able to 
compensate for the increased distortion energy. Whereas the TB model correctly predicts 
"k" to be the global minimum, the quantitative comparison with LSDA-GGA is rather poor. 
We suppose that this is a consequence of the small basis set of TB and more limited self 
consistency. Thus TB overestimates effects based on covalent bonding and lattice distortions 
because the Hamiltonian does not have the degrees of freedom of LSDA-GGA to find a lower 
energy in the case of an unfavorable structure. This is well illustrated in the case of "(110)" 
in figure 13(c). Clearly there is a large distortion of the bcc cubic unit cell surrounding the 
vacancy. The LSDA-GGA given the constraint of this distortion in the atomic structure is 
yet able to find a low energy electronic structure that can accommodate the constraint. The 
TB finds this much more difficult. 
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V. CONCLUSIONS 



The TB model presented here is intended as a physically better motivated and more 
transferable scheme as compared to the recently published orthogonal pd-mode\r^ or to 
existing classical interatomic potentials. Its transferability has been demonstrated using tests 
in both concentrated and dilute limits, for example successfully predicting the structure and 
energetics of cementite (section [HI Aj) and the migration path of C in 7-Fe (section [HI Bj) . It 
may be thought that this migration path particularly would expose the need for environment 
dependence in an empirical model^ 1 ^ Our model, using only two-center parameters, is able 
to deal with this through the use of an overlap matrix between non orthogonal Fe and C 
orbitals^ 

Our model also correctly describes the structure and energetics of the carbon dimer bound 
to a vacancy in a-Fe — a defect that is expected to take a central importance following the 
predictions of Forst et al— Apart from a large overestimation of the energy of the "(Oil)" 
dimer, our TB model properly orders the structures and predicts "k" to have the lowest 
energy although we were forced to modify an existing simple model for carbon in order to 
achieve the correct C2 bond length (see section Hi Bp . It is notable that published classical 
potentials^ and the minimal basis TB modeU^ cannot reproduce the stability of the carbon 
dimer. An exception is the recent classical potential of Lau et al— although this model 
greatly overestimates the binding energy of the (Oil) dimer. 

In view of the apparent significance of carbon dimers existing in the microstructure of 
steel and their possible interactions with hydrogen^ it is now a matter of importance that 
plausible and efficient quantum mechanical models are produced. From this point of view 
the present work assumes a particular significance. 
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